﻿using System;
using System.Collections.Generic;
using System.Linq;
using BSU.GraduateWork.Data.Model;
using BSU.Math.Algebra;
using BSU.Math.Geometry;

namespace BSU.GraduateWork.BL
{
    public class GridFreeMethod
    {
        #region Public properties

        public Problem Problem { get; set; }

        #endregion

        #region Private properties

        private List<Point4D> PrevU { get; set; }

        #endregion

        #region Constructors and Initializers

        public GridFreeMethod(Problem problem)
        {
            Problem = problem;
            Problem.Solution = new Solution
                                   {
                                       U = GenerateRegularPoints(Problem.InputParameters.N, Problem.InputParameters.X1, Problem.InputParameters.X2, Problem.InputParameters.Y1, Problem.InputParameters.Y2, Problem.InputParameters.Z1, Problem.InputParameters.Z2)
                                   };
        }

        public GridFreeMethod(int n,
                              double x1, double x2,
                              double y1, double y2,
                              double z1, double z2,
                              double g, double epsilon, double h, double alpha)
        {
            Problem = new Problem
                          {
                              InputParameters = new InputParameters
                                                   {
                                                       N = n,

                                                       X1 = x1,
                                                       X2 = x2,
                                                       Y1 = y1,
                                                       Y2 = y2,
                                                       Z1 = z1,
                                                       Z2 = z2,

                                                       G = g,
                                                       Epsilon = epsilon,
                                                       H = h,
                                                       Alpha = alpha
                                                   },
                              Solution = new Solution
                                             {
                                                 IterationsCount = 0
                                             }

                          };

            Problem.Solution.U = GenerateRegularPoints(Problem.InputParameters.N, Problem.InputParameters.X1, Problem.InputParameters.X2, Problem.InputParameters.Y1, Problem.InputParameters.Y2, Problem.InputParameters.Z1, Problem.InputParameters.Z2);
        }

        public static double F(double x, double y, double z)
        {
            return -2 * (y * (1 - y) * z * (1 - z) + x * (1 - x) * z * (1 - z) + x * (1 - x) * y * (1 - y));
        }

        #endregion

        #region Private methods

        private Matrix GetM(Point4D p, IList<Point4D> neighboringPoints)
        {
            int n = neighboringPoints.Count;

            var m = new Matrix(n + 1, 10);
            for (int i = 0; i < n; i++)
            {
                var point = neighboringPoints[i];

                m[i, 0] = 1;
                m[i, 1] = point.X - p.X;
                m[i, 2] = point.Y - p.Y;
                m[i, 3] = point.Z - p.Z;
                m[i, 4] = 0.5 * (point.X - p.X) * (point.X - p.X);
                m[i, 5] = (point.X - p.X) * (point.Y - p.Y);
                m[i, 6] = (point.X - p.X) * (point.Z - p.Z);
                m[i, 7] = 0.5 * (point.Y - p.Y) * (point.Y - p.Y);
                m[i, 8] = (point.Y - p.Y) * (point.Z - p.Z);
                m[i, 9] = 0.5 * (point.Z - p.Z) * (point.Z - p.Z);
            }
            m[n, 0] = 0; m[n, 1] = 0; m[n, 2] = 0; m[n, 3] = 0; m[n, 4] = 1;
            m[n, 5] = 0; m[n, 6] = 0; m[n, 7] = 1; m[n, 8] = 0; m[n, 9] = 1;

            return m;
        }

        private Matrix GetW(Point4D p, IList<Point4D> neighboringPoints)
        {
            int n = neighboringPoints.Count;
            var w = new Matrix(n + 1, n + 1);
            for (int i = 0; i < n; i++)
            {
                var point = neighboringPoints[i];
                double norm = (new Vector(point) - new Vector(p)).Norm();

                if (norm / Problem.InputParameters.H <= 1)
                {
                    w[i, i] = System.Math.Exp(-Problem.InputParameters.Alpha * norm * norm / (Problem.InputParameters.H * Problem.InputParameters.H));
                }
                else
                {
                    w[i, i] = 0;
                }
            }
            w[n, n] = 1;
            return w;
        }

        private Vector GetB(Point4D p, IList<Point4D> neighboringPoints)
        {
            int n = neighboringPoints.Count;
            var b = new Vector(n + 1);

            for (int i = 0; i < n; i++)
            {
                var point = neighboringPoints[i];
                b[i] = GetPoint(point, neighboringPoints).U;
            }
            b[n] = F(p.X, p.Y, p.Z);
            return b;
        }

        public void Solve()
        {
            var startTime = DateTime.Now;

            Problem.Solution.IterationsCount = 0;
            do
            {
                PrevU = new List<Point4D>(Problem.InputParameters.N);
                foreach (Point4D p in Problem.Solution.U)
                {
                    PrevU.Add(new Point4D(p));
                }

                for (int i = 0; i < Problem.Solution.U.Count; i++)
                {
                    Point4D point = Problem.Solution.U[i];
                    if (!point.IsBoundary)
                    {
                        var neighboringPoints = GetNeighboringPoints(point);

                        var m = GetM(point, neighboringPoints);
                        var w = GetW(point, neighboringPoints);
                        var b = GetB(point, neighboringPoints);

                        var mt = m.Transpose();
                        var tmp = mt * w;

                        var a = (tmp * m).Inverse() * tmp * b;
                        point.U = a[0];
                    }
                    else
                    {
                        point.U = Problem.InputParameters.G;
                    }
                }
                Problem.Solution.IterationsCount++;
            } while (!IsIterationsStopped());

            var endTime = DateTime.Now;

            Problem.Solution.ExecutionTime = endTime - startTime;
        }

        private bool IsIterationsStopped()
        {
            double summ1 = 0;
            double summ2 = 0;

            for (int i = 0; i < Problem.Solution.U.Count; i++)
            {
                summ1 += System.Math.Abs(Problem.Solution.U[i].U - PrevU[i].U);
                summ2 += System.Math.Abs(Problem.Solution.U[i].U);
            }
            return summ1 / summ2 < Problem.InputParameters.Epsilon;
        }

        private Point4D GetPoint(Point4D point, IEnumerable<Point4D> points)
        {
            var list = points
                .Where(p => p.X == point.X && p.Y == point.Y && p.Z == point.Z)
                .ToList();

            if (list.Count > 0)
            {
                return list[0];
            }
            return null;
        }

        private List<Point4D> GenerateRegularPoints(int count, double x1, double x2, double y1, double y2, double z1, double z2)
        {
            var points = new List<Point4D>(count);

            int boundaryPointsCount = (int)System.Math.Pow(count, 1.0 / 3);

            double hx = (x2 - x1) / boundaryPointsCount;
            double hy = (y2 - y1) / boundaryPointsCount;
            double hz = (z2 - z1) / boundaryPointsCount;

            for (int i = 0; i <= boundaryPointsCount; i++)
            {
                for (int j = 0; j <= boundaryPointsCount; j++)
                {
                    for (int k = 0; k <= boundaryPointsCount; k++)
                    {
                        points.Add(new Point4D(x1 + i * hx, y1 + j * hy, z1 + k * hz, 0,
                                               i == 0 || j == 0 || k == 0 ||
                                               i == boundaryPointsCount || j == boundaryPointsCount || k == boundaryPointsCount));
                    }
                }
            }

            return points;
        }

        private IList<Point4D> GetNeighboringPoints(Point4D point)
        {
            var list = Problem.Solution.U
                .Where(p => GetDistance(point, p) <= Problem.InputParameters.H &&
                            (p.X != point.X && p.Y != point.Y && p.Z != point.Z))
                .ToList();

            return list;
        }

        private static double GetDistance(Point4D a, Point4D b)
        {
            return System.Math.Sqrt((a.X - b.X) * (a.X - b.X) +
                                    (a.Y - b.Y) * (a.Y - b.Y) +
                                    (a.Z - b.Z) * (a.Z - b.Z));
        }

        #endregion

    }
}
